Energy absorption by "sparse" systems: beyond linear response theory 



(N 

O 

<N 



O 
(N 



X: 



Doron Cohen 

Departments of Physics, Ben-Gurion University of the Negev, P.O.B. 653, Beer-Sheva 84105, Israel 

The analysis of the response to driving in the case of weakly chaotic or weakly interacting systems 
should go beyond linear response theory. Due to the "sparsity" of the perturbation matrix, a resistor 
network picture of transitions between energy levels is essential. The Kubo formula is modified, 
replacing the "algebraic" average over the squared matrix elements by a "resistor network" average. 
Consequently the response becomes semi-linear rather than linear. Some novel results have been 
obtained in the context of two prototype problems: the heating rate of particles in Billiards with 
vibrating walls; and the Ohmic Joule conductance of mesoscopic rings driven by electromotive force. 
Respectively, the obtained results are contrasted with the "Wall formula" and the "Drude formula" . 
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I. INTRODUCTION 



This presentation concerns driven systems, like those 
illustrated in Figni whose dynamics is generated by an 
Hamiltonian that is represented by a matrix that has the 
generic structure 



n^^,,u = dmg{En} ~ f{t){Vnm} 



(1) 



Here i5„ are the ordered energy levels of the unperturbed 
system. Their density p[statcs/energy] is assumed to be 
roughly uniform. The system is driven by a low frequency 
stationary driving source f{t). The elements of the per- 
turbation matrix are Vnm- The induced transitions have 
rates that are proportional to |Kim P- We define 



X = {\Vn,nf} 



(2) 



Our interest concerns Hamiltonians that have a "sparse" 
perturbation matrix. This means that the majority of 
elements in X are small. To be more precise we assume 
that these elements have a log-wide distribution with a 
median that is much smaller compared with the aver- 
age. An example for such matrix is given in FigHJ and a 
typical histogram of the elements is presented in Fig ID 

The question that we ask is simple: Given X, what 
is the calculation that should be done in order to get 
the energy absorption rate (EAR). It makes sense that 
the result should be proportional to some weighted av- 
erage {{X)) over the matrix elements. Indeed, within 
the framework of linear response theory (LRT) the Kubo 
formula is doing just that - a weighted algebraic average. 

One should realize that the use of the Kubo formula for 
EAR calculation can be justified only in the very weak 
driving limit, provided there is a background "bath" that 
maintains quasi-equilibrium at any moment. However, if 
the driving is not very weak, compared with the relax- 
ation, then one should be worried: in order for the system 
to heat up, it is essential to have connected sequences of 
transitions, else the system is "stuck" . There is an obvi- 
ous analogy here with a resistor network calculation: due 
to the sparsity the energy absorption somewhat resem- 
bles a percolation process. 




FIG. 1: Model systems: a Billiard with a moving wall (up- 
per panel), and a Ring with a time dependent magnetic flux 
(lower panel). The deviation of the Billiard from intelligibility 
is quantified by a parameter u. It is due to the deformation 
of the boundary (as in the figure) or due to a deformation of 
the potential fioor (not illustrated). In the lower panel the 
Ring is regarded as a rectangular-like billiard with periodic 
boundary conditions in one of its coordinates. In the numer- 
ics the Ring has been modeled as a tight binding array of 
dimensions L x M. In the latter case the non-integrability 
was due to on-site disorder W . 



The bottom line of the above considerations is that the 
algebraic average of Kubo {{X))a, should be replaced by 
a resistor network average ((X))s, whose value is much 
smaller if the matrix is "sparse" . This should be regarded 
as an anomaly in the theory of response: it is an effect 
that arises upon quantization. Namely, one can charac- 
terize the sparsity of X by a parameter < s < 1 that is 
absent in the classical context, but has a dramatic effect 
in the quantum analysis. 
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FIG. 2: The upper panel is an image of the matrix X for 
a billiard very similar in shape to that of Fig[l] This matrix 
is "sparse". This can be deduced either by inspection, or 
by looking in the lower panel where the average value {x) of 
the elements (upper curves) is calculated as a function of the 
distance u) from the diagonal, and contrasted with the much 
smaller median values (lower curves) . For further details see 

[si,!!!. 

A few words are in order regarding the literature. We 
go here beyond the eonventional random matrix theory 
(RMT) perspective of ULB, because we are dealing 
with "sparse" matrices [JHI , possibly banded [7l-[l2j . 
Our view of LRT follows that of [H-ii. In the context 
of Billiards, LRT implies the "Wall formula" [Hl-i^, 
while in the context of meso scop ic conductance LRT 
implies the "Drude formula" [23] • In both cases one 
should take into account corrections that are related to 
correlations and level statistics. The quest for anomalies 
that cannot be explained by introducing corrections 
within the framework of LRT, but go beyond LRT, has 
some history [H, [H, IH-Ilil . The hue of study regarding 
the anomaly that arises due to "sparsity" is documented 
in (ssl - lioj . see acknowledgment. The resistor network 
analysis that is introduced below is inspired by [4l| - [46| . 
but generalizes its scope in a somewhat revolutionary 
way. 

Outline.- We first introduce with more details the 
model systems of FigQ] Then we outline the formalism 
of the EAR calculation, that is based on a simple Fermi 
golden rule (FGR) picture. Finally we present results 
that we have obtained for the dependence of the absorp- 
tion coefficient or the conductance on the sparsity, where 
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FIG. 3: Histogram of the in-band elements x of the matrix 
X for a rectangular billiard that has deformed potential floor. 
Different symbols refer to different values of the deformation 
parameter it. The vertical line is the average value of the 
elements, while the median is roughly at the locations of the 
peaks. The former unlike the latter is not sensitive to the 
degree of deformation. For details see [3^ . 

the latter is controlled by the degree of deformation or by 
the disorder in the system. Two appendices gives extra 
details on the resistor network calculation and on what 
we call "resistor network average" . 

II. THE MODEL SYSTEMS 

Assume that we have N non-interacting particles in a 
"box" . For presentation purpose assume a rectangular- 
like two dimensional billiard shaped box as in Fig.[T^, 
or optionally, imposing periodic boundary conditions in 
once direction, a ring shaped box as in Fig{T]3. The box 
is slightly deformed: either its walls are slightly curved, 
or optionally the potential floor is not flat, e.g. due to 
some scatterer or disorder. The system is driven by a 
low frequency stationary source. In FigU^ the driving is 
induced by moving a "piston" , while in Fig.[T]D it is by 
varying a magnetic flux through the ring. The Hamilto- 
nian matrix in the unperturbed energy basis takes that 
form of Eq.(IT|). In the case of the driven ring we have the 
identifications 

/(t) ^ m (3) 

(4) 

where $(i) is the magnetic flux, and Vnm are the matrix 
elements of the velocity operator. Note that by Faraday 
law — $ is the electromotive force (EMF). 

The driving induces transitions between energy levels. 
We assume stationary driving source, and define the its 
power spectrum as 

~S{lo) = FT (fit), f{0)) (5) 

Note that as far as EAR is concerned, it is a non-zero / 
that makes the Hamiltonian time dependent, correspond- 
ing to the EMF in the case of the ring. We assume low 
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frequency driving. Accordingly we write 



(6) 



where the prefactor is the RMS value that characterizes 
the driving intensity, and Sdc^) is a broadened delta func- 
tion whose line shape reflects the spectral content of the 
driving. 



III. THE ENERGY ABSORPTION RATE 

The driving induce transitions between energy levels, 
which implies diffusion in energy space. This diffusion is 
characterized by a coefficient [energy^ /time] for which 
we would like to have a formula. Assuming that D is 
known, the EAR is given by the following expression: 



EAR 



density x D 



(7) 



This is a straightforward generalization of Einstein type 
relations that are discussed in [l^ and in greater details 
in We can call it a diffusion-dissipation relation. 

What we label in Eq.© as "density" stands for the num- 
ber of particles (N) per energy, meaning N/T in the case 
of a Boltzmann occupation at tempeature T, or g in the 
case of a low temperature Fermi occupation. In the lat- 
ter case the role of the temperature is overtaken by the 
Fermi energy, namely g ^ N/E-p. 

What we would like to have, is a theory that allows 
the calculation of D. The formulas that we would like to 
advertise is 



D 



(8) 



Depending on the interpretation of ((|KimP)) this is the 
Kubo formula of LRT, or its resistor-network variation. 
In the latter case we refer to the results as the outcome 
of semi-linear response theory (SLRT). The latter term 
indicates that ((|T4mP))s unlike ((|KimP))a is a semi- 
linear rather than linear operation. The derivations of 
both the LRT and the SLRT variations of Eq. dH) are 
outlined in the next sections. 

In the case of an EMF driven Ring, it is convenient to 
re- write the EAR formulas Eq.© as 



EAR = G/2 



(9) 



where the so-called mesoscopic conductance is given by 
the expression 



G 



ng 



iilVn^nf)) (10) 



One can regard G as a mesoscopic version of the ab- 
sorption coefficient, while Eq.® can be regarded as the 
mesoscopic version of Joule law. 



IV. THE FERMI GOLDEN RULE PICTURE 
AND THE KUBO FORMULA 

The Hamiltonian in the standard basis is Eq.(IT]). We 
can transform it to the adiabatic basis: 



n = dia.g{En}- fit) 



En — E„ 



(11) 



The FGR transition rate from Em to En due to the low 
frequency noisy driving is: 



Vn 



En — En 



S {En — Em) 



(12) 



The FGR transitions lead to diffusion in energy space. 
Assuming that there is a background relaxation process 
that maintains at any moment quasi-equilibrium with oc- 
cupation probabilities p„ that are the same as in the ab- 
sence of driving, we get for the driving induced diffusion 



D 



n 

leading to the Kubo formula 

D = 



9 y^jEm-En)^ 



Cico)Sico)p 

ZTT 



(13) 



(14) 



with the spectral function 
Cicu) = FT {V{t)V{0)) 



(15) 



Y^P^Y. 27^^^^ " (Em-En)) 



Note that this spectral function reflects the band profile 
of X, as defined in Eq.Q, and illustrated in FigH 

It is important to realize that the Kubo formula 
Ea. (fT4l) is a linear functional of S{uj). The dependence 
on the matrix elements of X is linear too: 



D 



7r^p„|Km„P 5c{Em--En) 



(16) 



This can be formally written as Eq. with an implied 
definition of the weighted algebraic average ((|V^mP))a- 
For completeness we note that in the context of con- 
ductance calculation the popular textbook version of the 
Kubo formula is 



St(E„-Ef) SAE,n-E„) (17) 



71, m 



This expression is implied by Ea. (|16p , with averaging over 
the levels in the vicinity of the Fermi energy. Namely 
Pn ~ Q'^^riEn—Ep), where the width of S^O is deter- 
mined by the temperature. 
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FIG. 4: The driving induces transitions between levels E„ 
of a closed system, leading to diffusion in energy space, and 
hence an associated heating. The diffusion coefficient D can 
be calculated using a resistor network analogy. Connected 
sequences of transitions are essential in order to have a non- 
vanishing result, as in the theory of percolation. 



V. THE CALCULATION OF THE DIFFUSION 
COEFFICIENT 

We would like to consider circumstances in which the 
driving induced transitions are faster compared with the 
background relaxation. In such circumstances the occu- 
pation probabilities p„ are no longer as in equilibrium. 
For the purpose of analysis we simply neglect the bath, 
and describe the dynamics by a rate equation 



—T- = - 2^ Wnm{Pn " Pmj 



(18) 



where the rates Wnm sue determined by Eq. (1121) . The 
matrix w = {wnm} can be regarded as a quasi one- 
dimensional network. See Fig.lH Optionally one may 
interpret the rate equation as a probabilistic description 
of a random walk process, where Wnm is the probability 
to hop from to to n per unit time. The local spreading 
is described by 



Var(n) 



^[w„,not] (n - no) 



(19) 
(20) 



It follows that the course-grained spreading should be 
described by a diffusion equations 



dpn _ jj d'^ 

dt d-n? 



(21) 



where n is regarded as a continuous variable. The diffu- 
sion equation is formally a continuity equation 

% - -^In (22) 

ot on 
where the current is given by Pick's law: 



'D-^p„ 
on 



(23) 



If we have a sample of length N then 

I = X[PN -Po] 



(24) 



This shows that D/N is formally like the inverse resis- 
tance of the chain: it is the ratio between the current and 
the "potential difference" . We therefore can use standard 
recipes of electrical engineering in order to calculate its 
value. For example, if we have only near-neighbor tran- 
sitions then "adding connectors in series" implies 



D 

N 



■ N 

E: 

.n=l 



1 



(25) 



In general we use the notation D — [[w]], where the 
doubled brackets stand for inverse resistivity calculation, 
as discussed in AppEl 

In the above analysis we have assumed unit distance 
between sites. If the mean level spacing is g~^, the ex- 
pression for the diffusion coefHcient should be re-scaled 
as follows: 



Recall that the 



D = g- 



D = g''[[w]] (26) 
im are given by Eq. ([T2l) . it follows that 



En ~ En 



2Tr\f\^Sc{E„-E„,) 



(27) 



leading to Eg. pH)) , with an implied definition of the resis- 
tor network average ((|KmiP))s- The definition and the 
calculation of the latter are further discussed in ApplBl 



VI. THE WALL FORMULA AND BEYOND 

The roughest estimate for the diffusion that is induced 
by a vibrating wall is known as the "Wall formula" |21^ 
l2a |. Its original derivation is based on a simple kinetic 
picture: it is based on the assumption that collisions with 
the vibrating wall are not correlated. This leads in the 
two dimensional case [24, 50l to the result 



Do 



P 



(28) 



where m is the mass of the particles, and is the linear 
dimension of the box as illustrated in Fig{T] The result 
assumes a microcanonical preparation at energy E, and 
we have defined v-e ~ \2E / m]^/"^ . If we have a Boltz- 
mann occupation, the expression should be averaged ac- 
cordingly. If we have low temperature Fermi occupation, 
what counts in Eq.© is the the value at E ~ Ep. 

Within the framework of LRT the same result is ob- 
tained from the Kubo formula Eg. provided C{Ld) 
is flat, i.e. provided there are no correlations between 
collisions. In practice there are correlations leading to 
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D = gcDo with gc that can be either smaller or larger 
than unity depending on the geometry. In the quan- 
tum calculation gc is slightly affected by the level spacing 
statistics. 

Within the framework of SLRT one has to calculate the 
resistor network average of the X matrix. If this matrix 
is sparse, the result becomes much suppressed, leading to 
D = gsgcDo, where 



gs 



{{\Vnm\'))s 
((|K,„P))a 



(29) 



We emphasize that g^ reflects an anomaly: it depends 
on the sparsity s of the matrix, a parameter that has 
no meaning in the classical context. For more details, 
including RMT analysis of this dependence see [s^ ■ 

Some numerical results for gc (LRT) and g = gsgc 
(SLRT) arc presented in FiglH The calculation is done 
with the C{uj) of FigJ21 The h dependence of the LRT 
"classical" result is due to the lower cutoff of the doj in- 
tegral in the Kubo formula Eq. (fT4|) . that is sensitive to 
mean level spacing. If C{uj) were "flat" the result would 
not be much sensitive to ft. As expected the SLRT cal- 
culation gives a much smaller result. The effect becomes 
more conspicuous for smaller deformations, for which the 
sparsity is "stronger" . 
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VII. THE DRUDE FORMULA AND BEYOND 

The EAR by particles that are driven by an oscillat- 
ing electric field, due to an induced EMF, is very similar 
problem to that of particles driven by an oscillating wall. 
Here the simple result that is based on kinetic considera- 
tions is known as the "Drude formula" . As in the case of 
the "Wall formula" it is assumed that scattering events 
are uncorrclated, leading to the estimate 



(30) 



where e is the charge of the electron, L is the length of 
the ring, and £ is the mean free path between collisions. 
Note that T> = v^i is the spatial diffusion coefficient, 
while (e<f>/L) is the energy that is gained per circulation. 
The implied result for the conductance can be written as 



Go 



M - 

27rft L 



(31) 



where M = mv-ELy/-K is the number of open modes. This 
way of writing the Drude formula is very illuminating 
because it refiects Ohm law, and the units are as in the 
Landauer formula. For clarity we have restored the ft in 
this formula. 

Within the framework of LRT, the same result is ob- 
tained from the Kubo formula Ea.([T4)). provided C{uj) of 
Ea. (|15p is a Lorentzian that reflects an exponential decay 
of the velcoity-vclocity correlation function. In practice 
there are extra correlations leading to D — gcDo with gc 



FIG. 5: The scaled absorption coefficient gc (LRT) and 
g = QsQc (SLRT) versus the dimensionless 1/h (upper panel), 
and versus the dimensionless deformation parameter u (lower 
panel). The input for this analysis is the matrix X of Fig[2] 
The calculation of each point has been carried out on a 
100 X 100 sub-matrix of X centered around the fi implied 
energy E. The "untextured" data points are calculated for 
an artificial random matrices with the same bandprofile and 
sparsity. The complementary lower panel is oriented to show 
the small u dependence within an energy window that corre- 
sponds to 1/fi ~ 9. For further details see [39l.l40l|. 



that can be either smaller or larger than unity depending 
on the geometry. In the quantum calculation gc is slightly 
affected by the level spacing statistics, and the correction 
is of order (guic)~^- This is sometimes regarded as a vari- 
ation of weak localization corrections |27| . 

Within the framework of SLRT one has to calculate the 
resistor network average of the matrix { | P } ■ Here one 
should distinguish between different regimes, depending 
on the strength W of the disorder. In the Born ap- 
proximation the mean free path is £ cx W'^ , while the 
localization length is £oo ~ Mi- The diffusive regime, 
where there is no issue of sparsity, requires an interme- 
diate strength of disorder, such that £ ^ L ^ £oo- For 
either stronger or weaker disorder, the matrix {|w,i,n|^} 
becomes "sparse". This is because the eigenstates be- 
come non-ergodic: either they are localized in real space 
(for strong disorder) or in mode space (for weak disor- 
der). Note also that for very weak disorder ("clean" ring) 
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FIG. 6: Schematic illustration that sketch the dependence 
of the DC mesoscopic conductance on the strength of the 
disorder. It should be regarded as a caricature of FiglT] The 
level width V = fkOc affects the so-called weak localization 
correction in the diffusive regime. In the other two regimes 
of either weak or strong disorder, the perturbation matrix 
becomes "sparse" and consequently G is suppressed compared 
with Drude. 



each cigcnstate occupies a single mode (up to small cor- 
rection). For detailed analysis that supports the above 
picture see [svj . 

The expectation with regard o the dependence of G/Gq 
on the strength of the disorder is sketched in Fig[6l Some 
numerical results for both the LRT and the SLRT con- 
ductance are presented in FiglT] The calculation is done 
for a tight binding model. The SLRT result in the An- 
derson localization regime is completely analogous to the 
reasoning of variable range hopping [4ll445j , as explained 
in [4^ . It should be appreciated that in our approach all 
regimes are treated on equal footing. 

In the ballistic regime, contrary to the Drude expecta- 
tion, the conductance becomes worse as the disorder is 
reduced. This looks strange, but can easily be rational- 
ized if we think about the extreme case of no disorder: in 
the absence of scattering the particle stays all the time 
in the same mode; hence irreversible diffusive spreading 
in energy is impossible. 



VIII. CONCLUSIONS 

The random matrix approach of Wigner ('^ 1955) is 
based on the observation that in generic circumstances 
the perturbation can be represented by a random matrix 
whose elements are taken from a Gaussian distribution. 
Our interest in this presentation concerned a restricted 
class of "sparse" systems for which this observation does 
not hold. In such weak quantum chaos circumstances 
the elements are characterized by a log- wide distribution. 
Consequently, the response, and in particular the energy 
absorption, are similar to a percolation process, and their 
analysis requires a novel resistor-network approach. 

Besides the quantitative issue, the experimental finger- 
print of the resistor-network calculation is the implied 
semi-linearity of the response. In the SLRT regime, i.e. 
if the driving is the predominant mechanism for transi- 
tions between levels, one expects the combined effect of 
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FIG. 7: Various measures of sparsity (upper panel), and 
the mesoscopic conductance (lower panels) as a function of 
the disorder W in the Ring. The vertical line separate be- 
tween the clean, ballistic, diffusive and localization regimes. 
Note that the scaled conductance in arbitrary units equals 
((bmnp>>. The Drude, the LRT and the SLRT results are 
displayed in both normal and log scale. We see that in the 
ballistic regime the SLRT conductance becomes worse as the 
disorder becomes weaker, in opposition with the Drude ex- 
pectation. For further details see [3^ . [stI ] 

two independent sources to be super-linear, namely, 

D[Sa{L0) + Sb{uj)\ > D[Sa{i0)] + D[Sb{uj)] (32) 

but still scmi-hnear D [A^(a;)] ^ X D [S{uj)] . We have 
provides in this presentation two prototype examples 
where an SLRT anomaly can arise: heating of particles 
that are trapped in billiards with vibrating walls; and 
Joule heating of charged carriers that are driven by an 
induced electro-motive force. 
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Appendix A: The resistor network calculation 

In this appendix we explain how the inverse resis- 
tivity G = [[G]] of a one-dimensional resistor network 
G = {Gnm} is calculated. We use the language of electri- 
cal engineering for this purpose. In general this relation 
is semi-liner rather than linear, namely [[AG]] = A[[G]], 
hnt[[A + B]]^[[A]] + [[B]]. 

There are a few cases where an analytical expression is 
available. If only near neighbor nodes are connected, al- 
lowing Gn,n+i = On to bc different from each other, then 
"addition in series" implies that the inverse resistivity 
calculated for a chain of length N is 



G 



(Al) 



If Gnm ~ gn~m IS a function of the distance between the 
nodes n and m then it is a nice exercise to prove that 
"addition in parallel" implies 



G 



(A2) 



In general an analytical formula for G is not available, 
and we have to apply a numerical procedure. For this 
purpose we imagine that each node n is connected to a 
current source /„ . The Kirchhoff equations for the volt- 
ages are 



(A3) 



This set of equation can be written in a matrix form: 

GV = I (A4) 

where the so-called discrete Laplacian matrix of the net- 
work is defined as 



G„ 



G„ 



(A5) 



This matrix has an eigenvalue zero which is associated 
with a uniform voltage eigenvector. Therefore, it has 
a pseudo-inverse rather than an inverse, and the Kirch- 
hoff equation has a solution if and only if X^n ^" ~ '^^ 
order to find the resistance between nodes n^^ = and 
7i„„, = N, we set Iq ~ 1 and /jv = — 1 and /„ = other- 
wise, and solve for Vq and Vn- The inverse resistivity is 
G=[iVo-VN)/N]-\ 



Appendix B: The resistor- newtworli average 

We use the notation {{X)) in order to indicate the 
weighted average value of its elements. First we would 
like to define the standard algebraic average. It is essen- 
tial to introduce a weight function that defines the band 



of interest. In the physical context this function reflects 
the spectral content of the driving source. Namely, we 
define F{r) as the normalized version of S{u>), such that 
X^r ^('') ~ 1' where r = n — m is the energy difference 
uj = En — Em in integer units. The bandwidth in these 
dimensionless units (be = qloc) it is assumed to be quan- 
tum mechanically large {be 3> 1). The algebraic average 
is defined in the standard way: 



((^))a 



-y 



F{n — m) Xn 



(Bl) 



where N is the size of the matrix, which is assumed to be 
very large. The algebraic average is a linear operation, 
meaning that 

((AX)) = A((X)) (B2) 
{{X + Y)) = {{X)) + {{Y)) (B3) 

There arc different type of "averages" in the literature, 
such as the harmonic average, geometric average, and we 
can also include the median in the same list. All these 
"averages" are semi-linear operations because only the 
((AX)) = A((X)) property is satisfied for them. Irre- 
spective of the semi-linearity issue any type of average 
should satisfy the following requirement: if all the ele- 
ments equal to the same number, then also the average 
should equal the same number. 

In this presentation we highlight a new type of average 
that we call a resistor-network average: 



2F{n-m) 



Xn 



{n — my 



(B4) 



Writing the expression above as [[wnm]], one should re- 
alize that the Wnm can be regarded as FGR transition 
rates. Using Ea. (|A2p it is not difficult to show that if all 
the elements X„,„ are the same number, then also their 
resistor network average is the same number. While in 
general 



{{X))s < {{X))a 



(B5) 



Typically the resistor network average is bounded from 
below by the median. In order to get a realistic estimate 
in the case of a "sparse" matrix one can use a generalized 
variable range hopping scheme that we have developed 
in [13. 
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